Why use R for neuroimaging?

R has excellent tools for development and reproducibility

Development and maintenance is independent of “elite” and “esoteric” user groups.
- unlike pretty much all other neuroimaging software tools - R plays nicely with Shell scripts, Python, C++, etc.

RStudio and RMarkdown

  • We can easily create new functions, scripts, and dynamic documents that specify and describe our neuroimaging workflow.
  • RStudio is a feature-rich IDE with everything you need to implemen a workflow or develop your own.
  • With RMarkdown, we can organize the descriptsions of your workflow, all of your code (including, R, bash, python, C++, etc.) in one place, AND we can describe our results and make publication ready figures and tables in a convenvient, shareable PDF or HTML format.

Installation

Linux makes the world a more useful place

Many neuroimaging tools require some flavour of Linux.

FSL and ANTs require Linux (and so does R if you want to use these tools).

Remember:

Neuroimaging data is not special, its just numbers!

It’s just a lot of numbers…

Numbers in specific locations.

When dealing with neuroimaging data in R, we encounter a few important issues with BIG Data

  • R has limits to the number of values that an object can contain
.Machine$integer.max
## [1] 2147483647
  • Big data can require BIG RAM
  • R likes to hold on to the RAM its been given, even after it should have let it go.
  • We can force R to give RAM back to the operating system when we want
# It is good practice to clear your workspace when starting a new project or workflow
rm(list=ls()) # clear all objects from the R environment
gc()# force R to release RAM to the OS
##          used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 487492 26.1     940480 50.3   940480 50.3
## Vcells 920922  7.1    1650153 12.6  1112749  8.5

R is extended by adding PACKAGES

Adding packages to your R installation is easy!

install.packages("ggplot2")
pkgs <- c("devtools", "lme4", "lmerTest", "fslr")
install.packages(pkgs)
# Install ANTsR
devtools::install_github("stnava/cmaker")
devtools::install_github("stnava/ITKR")
devtools::install_github("stnava/ANTsR")

# Install TKoscik's Tools
devtools::install_github("TKoscik/tk.nifti") # name change pending / CRAN
devtools::install_github("TKoscik/tk.r.misc") # name change pending / CRAN
devtools::install_github("TKoscik/cluster.nii") # coming soon
devtools::install_github("TKoscik/mRi") # TBD
source("/path/to/your/function.R")

Datatypes and Objects

In R your data can take a few different forms

Data are stored in the R environment as Objects

Vector: 1-dimensional object of the same datatype
Matrix: 2-dimensional object of the same datatype
Array: N-dimensional object of the same datatype

R has some more powerful datatypes

Dataframe: 2-dimensional object that allows mixed datatypes
List: 1-dimensional object that contains other objects

Data I/O

R has a suite of utilities to bring your data in

  • This includes Excel spreadsheets, as well as SPSS, SAS, and Matlab files
  • reading data from binary files, bit-by-bit or line-by-line is also relatively simple
  • unless its freesurfer data, which uses 3-bit integers in their header

Delimited files simplify data storage

  • comma- or tab-delimited data

use the function read.csv()

Indexing

Square brackets, [ ] are used to index your objects

a <- sample(1:10, 5, replace = TRUE)
print(a)
## [1] 8 7 1 2 6
a[3]
## [1] 1
a[c(1,3,5)]
## [1] 8 1 6
df <- data.frame(a = 1:5, b = rnorm(5))
df[2, ]
##   a         b
## 2 2 0.6388074
df[1,"b"]
## [1] 1.065538
df$a[1:3]
## [1] 1 2 3
temp.ls <- list(a=a, df=df)
temp.ls[[1]]
## [1] 8 7 1 2 6
temp.ls$a
## [1] 8 7 1 2 6

Operators

Operators are the symbols that indicate the operation to be implemented

Assignment Operators

Operator Operation
<- leftward assignment to current environment
= leftward assignment (within current context only)
-> rightward assignment
<<- or ->> parent environment assignment

<- and = are almost interchangeable, however using <- when assigning into a variable the element is sent to the environment as a standalone object as well

Arithmetic operators

Operator Operation
+ addition
- subtraction
* multiplication
/ division
^ or ** exponentiation
%% modulus
%/% integer division

Logical Operators

Operator Operation
== exactly equal to
< less than
> greater than
<= less than or equal to
>= greater than or equal to
!= not equal to
! not
|| OR
| Element-wise OR
&& AND
& Element-wise AND

Functions

To do things in R you call functions
* use parentheses ( ) to call a function
* assign output of a function to a variable with <-
* functions have arguments to control what they do

*** Example HERE

Modelling Syntax

One of the most useful concepts in R that is the statistical formula notation that most functions use

Basic Modelling syntax

dependent variables are modelled as a function of independent variables, with a ~

DV ~ IV1 terms: IV1

Main effects are added to the model with a +

DV ~ IV1 + IV2 terms: IV1, IV2

Interaction terms are added with :

DV ~ IV1 + IV2 + IV1:IV2 terms: IV1, IV2, IV1:IV2

Main effects and interactions can be added simultaneously with *

DV ~ IV1 * IV2 terms: IV1, IV2, IV1:IV2

higher-order terms are added by combining, multiple : or ’*’

DV ~ IV1 * IV2 * IV3 terms: IV1, IV2, IV1:IV2, IV1:IV3, IV2:IV3, IV1:IV2:IV3

parentheses provide control over terms in your model

DV ~ (IV1 + IV2) * IV3 terms: IV1, IV2, IV1:IV3, IV2:IV3

Practical

  1. structural preprocessing
    • bias correction
    • brain extraction
    • registration to template
  2. functional preprocessing
    • motion correction
    • registration to structural
    • data cleaning
    • temporal filtering
    • spatial smoothing
  3. modelling
    • task
    • resting-state
  4. clusterizing
  5. visualization

1. Structural preprocessing

Set Environment

rm(list=ls())
gc()

# Load ANTs library
library("ANTsR")

Bias correction - N4 method

Perform bias correction before brain extraction to improve extraction performance. We will redo this using our brain mask to initilize it to get better results after

data.dir <- paste0(getwd(), "/data/debias") # Directory for raw structural data
save.dir <- paste0(getwd(), "/data/debias") # Directory to save processed data

# load raw anatomical image
img.raw <- antsImageRead(filename = paste0(data.dir, "/T1.nii.gz"))

# preform N4 bias correction
img.n4 <- abpN4(img = img.raw,
                intensityTruncation = c(0.025, 0.975, 256),
                mask = NA,
                usen3 = FALSE)

# save N4 corrected image
antsImageWrite(image = img.n4,
               filename = paste0(save.dir, "/T1_n4.nii.gz"))

Brain extraction

template.dir <- paste0(getwd(), "/data/templates") # Directory for registration and brain extraction templates

# Load brain extraction template
tem <- antsImageRead(paste0(template.dir, "/T_template0.nii.gz"))
# Load registration mask
temmask <- antsImageRead(paste0(template.dir, "/T_template0_BrainCerebellumProbabilityMask.nii.gz"))

# perform brain extraction
img.bex <- abpBrainExtraction(img = img.n4,
                              tem = tem,
                              temmask = temmask)

# Save extracted brain
antsImageWrite(image = img.bex[[1]],
               filename = paste0(save.dir, "/T1_bex.nii.gz"))
antsImageWrite(image = img.bex[[2]],
               filename = paste0(save.dir, "/T1_bmask.nii.gz"))

Redo bias correction

# preform N4 bias correction
img.n4.bex <- abpN4(img = img.bex[[1]],
                intensityTruncation = c(0.025, 0.975, 256),
                mask = NA,
                usen3 = FALSE)

# save N4 corrected image
antsImageWrite(image = img.n4,
               filename = paste0(save.dir, "/T1_n4_bex.nii.gz"))

Register structural to common space template

# Load template image
fixed.img <- antsImageRead(paste0(template.dir, "/MNI152_T1_2mm_brain.nii.gz"))

# Load brain extracted T1
moving.img <- img.n4.bex

# perform rigid, linear, and diffeomorphic (nonlinear) warping
tform <- antsRegistration(fixed = fixed.img,
                            moving = moving.img,
                            typeofTransform = "SyN",
                            initialTransform = NA,
                            outprefix = paste0(save.dir, "/t1_tform_"))

# apply warping
img.reg <- antsApplyTransforms(fixed = fixed.img,
                               moving = moving.img,
                               transformlist = tform$fwdtransforms,
                               interpolator = "linear",
                               imagetype = 2)

# save registered image
antsImageWrite(image = img.reg, filename = paste0(save.dir, "/T1_reg.nii.gz"))

2. Functional preprocessing

Motion correction

func.dir <- paste0(getwd(), "/data/func") # specify directory for raw functional data

func.img <- antsImageRead(paste0(func.dir, "BOLD_run1.nii.gz")) # specify functional data
func.mean <- getAverageOfTimeSeries(func.img)

mocor <- antsMotionCalculation(img = func.img,
                               fixed = func.mean,
                               moreaccurate = TRUE,
                               framewise = TRUE)

mocor.params <- cbind(mocor$moco_params[ ,3:ncol(mocor$moco_params)], mocor$fd)

antsImageWrite(image = mocor$moco_img,
               filename = paste0(save.dir, "/BOLD_run1_mocor.nii.gz"))
write.table(x = mocor.params,
            file = paste0(save.dir, "/BOLD_run1_mocor.csv"),
            quote = FALSE, row.names = FALSE, col.names = TRUE, sep = ",")

Registration to structural

Register mean BOLD image to structural image

# load fixed image, and resample to 2mm isotropic space

fixed.img <- resampleImage(antsImageRead(paste0(save.dir, "/T1_n4.nii.gz")), rep(2,3))
moving.img <- getAverageOfTimeSeries(paste0(save.dir, "/BOLD_run1_mocor.nii.gz"))
        
bold.tform <- antsRegistration(
  fixed = fixed.img,
  moving = moving.img,
  typeofTransform = "SyNBold",
  initialTransform = NA,
  outprefix = paste0(save.dir, "bold_tform_"))

bold.to.t1 <- antsApplyTransforms(
          fixed = fixed.img,
          moving = moving.img,
          transformlist = bold.tform$fwdtransforms,
          interpolator = "linear",
          imagetype = 3)

Register to template

bold.to.template <- antsApplyTransforms(
  fixed = antsImageRead(paste0(template.dir, "/MNI152_T1_2mm_brain.nii.gz")),
  moving = bold.to.t1,
  transformlist = tform$fwdtransforms,
  interpolator = "linear",
  imagetype = 3)

antsImageWrite(image=bold.to.template, filename=paste0(save.dir, "BOLD_run1_norm.nii.gz"))

Data cleaning

img.bold <- antsImageRead(paste0(save.dir, "BOLD_run1_norm.nii.gz"))
img.bold.mean <- getAverageOfTimeSeries(img.bold)
img.bold.mask <- getMask(img.bold.mean)

bold.mx <- timeseries2matrix(img.bold, img.bold.mask)
dvars.pre <- computeDVARS(bold.mx)
gc()
compcor.vars <- compcor(fmri = img.bold,
                        mask = img.bold.mask,
                        ncompcor = 6,
                        variance_extreme = 0.975,
                        randomSamples = 1,
                        returnv = FALSE,
                        returnhighvarmat = FALSE,
                        returnhighvarmatinds = FALSE,
                        highvarmatinds = NA,
                        scale = FALSE)
gc()

bold.mx <- residuals(lm(bold.mx ~ scale(compcor.vars)))
gc()

img.bold <- matrix2timeseries(referenceImage = img.bold,
                              maskImage = img.bold.mask,
                              timeSeriesMatrix = bold.mx)

antsImageWrite(image = img.bold,
               filename = paste0(save.dir, "/BOLD_run1_compcor.nii.gz"))

write.table(x = compcor.vars,
            file = paste0(save.dir, "/BOLD_run1_compcor.csv"),
            quote = FALSE, row.names = FALSE, col.names = TRUE, sep = ",")

Temporal filtering

img.bold <- antsImageRead(paste0(save.dir, "/BOLD_run1_compcor.nii.gz"))
img.bold.mean <- getAverageOfTimeSeries(img.bold)
img.bold.mask <- getMask(img.bold.mean)
bold.mx <- timeseries2matrix(img.bold, img.bold.mask)

bold.mx <- frequencyFilterfMRI(boldmat = bold.mx,
                               tr = antsGetSpacing(img.bold)[4],
                               freqLo = opts$tfilter.freqLo,
                               freqHi = opts$tfilter.freqHi,
                               opt = opts$tfilter.opt)

img.bold <- matrix2timeseries(referenceImage = img.bold,
                              maskImage = img.bold.mask,
                              timeSeriesMatrix = bold.mx)

antsImageWrite(image = img.bold,
               filename = paste0(save.dir, "/BOLD_run1_tfilter.nii.gz"))

antsImageWrite(image = img.bold,
               filename = paste0(save.dir, "/BOLD_run1_prep.nii.gz"))

3. Modelling - Task fMRI

Deconvolving BOLD - Single-trial Betas

The temporal sampling of the BOLD signal (one volume at each TR) is usually different from the onsets of your task trials Thus, we want to get the BOLD signal into the space of our task trials. We do this by “Deconvolving” the BOLD signal, which is estimating the strength of the BOLD signal for each task trial.

Canonical HRF

precision <- 2
params <- list()
params$precision <- 10^precision
params$kernel.length <- 32
params$a1 <- 0.0083333333
params$a2 <- 1.274527e-13
params$p1 <- 5
params$p2 <- 15
params$t <- seq(from = 0, to = params$kernel.length, by = 1/params$precision)
params$hrf <- exp(-params$t) * (params$a1 * params$t^params$p1 - 
                                  params$a2 * params$t^params$p2)
params$hrf <- params$hrf[3:length(params$hrf)]
params$hrf <- params$hrf/max(params$hrf)
    
plot(params$hrf)

Deconvolve timeseries to put BOLD signals into the space of task data

ts <- matrix(0,
             nrow=(n.tr * tr + params$kernel.length) * params$precision,
             ncol = nrow(tf))
for (i in 1:nrow(tf)) {
  ts[round(tf$onset[i] * params$precision + 1),i] <- tf$ev[i]
  ts[ ,i] <- convolve(ts[ ,i], rev(params$hrf), type = "open")[1:nrow(ts)]
}
conv.ts <- ts[seq(from = 1, to = (n.tr * tr * params$precision), by = tr * params$precision), ]

uf <- data.frame(sim.fmri = rnorm(n.tr, 0,1)*25,
                 conv.ts)
uf.melt <- melt(uf[ ,c(1,10,20,30)])
uf.melt <- data.frame(1:149, uf.melt)

ggplot(uf.melt, aes(x=X1.149, y=value, color=variable)) +
  theme_bw() +
  geom_line(size=1)

deconv.mdl <- lm(sim.fmri ~ ., uf)
tf$stbs <- coef(deconv.mdl)[2:(nrow(tf)+1)]

plotf <- melt(tf[ ,c(7,16)])
plotf <- data.frame(trial = 1:36, plotf)
ggplot(plotf, aes(x=trial, y=value, color=variable)) + theme_bw() +
  geom_line(size=1)

Modelling

rm(list=ls())
gc()
library("mRi")
library("parallel")

data.4d <- "/full/path/to/directory/containing/4d/files/or/a/list/of/file/names"
data.paradigm <- "/full/path/to/paradgm.csv"
data.mask <- "/full/path/to/regoin/of/interest/mask(s).nii"
save.dir <- "/empty/directory/to/save/results"
file.prefix <- "some.name"
num.cores <- detectCores()-1

my.command <- "my.function(df, coords, img.dims, save.dir, prefix)"
my.function <- function(df, coords, img.dims, save.dir, prefix) {
  mdl1 <- lm(fmri ~ ev, df)
  mdl1.coef <- as.data.frame(summary(mdl1)$coef)
  table.to.nii(mdl1.coef, coords, img.dims, save.dir, prefix, mdl1)
  
  df$choice <- (df$response + 1) / 2
  mdl2 <- glm(choice ~ fmri * ev, df)
  mdl2.coef <- as.data.frame(summary(mdl2)$coef)
  table.to.nii(mdl2.coef, coords, img.dims, save.dir, prefix, mdl2)
}
run.voxel(data.4d, data.pradigm, my.function, my.command, save.dir, file.prefix, num.cores, vx.order = "rand", verbose = TRUE)
  1. Modelling - Resting-State fMRI
rm(list=ls())
gc()
library("mRi")
library("parallel")

data.4d <- "/full/path/to/directory/containing/4d/files/or/a/list/of/file/names"
seed.vxl <- c(40,40,40)
data.paradigm <- data.frame(seed = load.voxel(parse.4d(data.4d)[[1]], coords = seed.vxl))
data.mask <- "/full/path/to/regoin/of/interest/mask(s).nii"
save.dir <- "/empty/directory/to/save/results"
file.prefix <- "some.name"
num.cores <- detectCores()-1

my.command <- "my.function(df, coords, img.dims, save.dir, prefix)"
my.function <- function(df, coords, img.dims, save.dir, prefix) {
  temp <- cor.test(df$fmri, df$seed)
  cor.result <- data.frame(r = temp$estimate,
                           p = temp$p.value,
                           t = temp$statistic,
                           df = temp$parameter)
  table.to.nii(cor.result, coords, img.dims, save.dir, prefix)
}
run.voxel(data.4d, data.pradigm, my.function, my.command, save.dir, file.prefix, num.cores, vx.order = "rand", verbose = TRUE)

4. Clustering

rm(list=ls())
gc()

library(mRi)
cluster.valley(tmap = "",
               tvol = 2,
               pmap = "",
               pvol = 2,
               mask = "",
               alpha = 0.05,
               min.size = 10,
               save.dir = "", 
               file.prefix = NULL)

5. Visualization

R has some fantastic visualization tools to make figures ready for publication

rm(list=ls())
gc()

library("mRi")
library("tk.r.misc")
library("ggplot2")
library("gridExtra")
library("raster")

overlay.plot <- draw.overlay(anat.nii = "/full/path/to/anatamoical/image/for/background.nii",
                             over.nii = "/full/path/to/nii/with/values/to/overlay.nii",
                             over.vol = 2,
                             over.color = list(good.color("cool")),
                             mask.nii = "/full/path/to/mask.nii",
                             mask.vol = 1,
                             roi.nii = "/full/path/to/roi/map.nii",
                             roi.val = 1:10,
                             roi.color = "#ff64ff",
                             orientation = "sagittal",
                             save.plot = FALSE,
                             return.plot = TRUE)
grid.arrange(overlay.plot, newpage = TRUE)